*************************************************
* Purpose -- TWFE analysis using QCEW data
************************************************
clear all
* RUN 00_path_master.do FIRST TO GET FILEPATHS


************************************************
* Loading different samples
************************************************
* 1. CZ-state
frame create twfe1
frame change twfe1

use "$clean/qcew/cz_state_sample.dta", clear

* creating frames for sub-periods

frame copy twfe1 twfe1_00
frame copy twfe1 twfe1_16

frame twfe1_00: keep if inrange(year,2000,2016)
frame twfe1_16: keep if year <= 2016


* 2. MSCZ

frame create mscz
frame change mscz

use "$clean/qcew/mscz_sample.dta", clear

* creating frames for sub-periods
frame copy mscz mscz_00
frame copy mscz mscz_16

frame mscz_00: keep if inrange(year,2000,2016)
frame mscz_16: keep if year <= 2016


* 3. BCP

frame create bcp
frame change bcp

use "$clean/qcew/bcp_sample.dta", clear

* creating frames for sub-periods
frame copy bcp bcp_00
frame copy bcp bcp_16

frame bcp_00: keep if inrange(year,2000,2016)
frame bcp_16: keep if year <= 2016


************************************************
* Creating different frames for DL regressions
************************************************

* getting MW data for 1985-1989 and 2020-2021 for 1990-2019 time period lags and leads

frame copy bcp bcp_dl
frame change bcp_dl
preserve
	use "$raw/other/mw_state_annual_vz.dta", clear
	gen mw = max(max_mw, max_fed_mw)
	keep statefips year mw
	keep if inrange(year,1985,1989) | inrange(year,2020,2022)
	
	tempfile mw
	save `mw'
restore
 

forvalues i = 1985/1989 {
	insobs 1
	replace year = `i' if _n == _N
}

forvalues i = 2020/2022 {
	insobs 1
	replace year = `i' if _n == _N
}

fillin pair_id_county year
drop if mi(pair_id_county)
bys pair_id_county (year): replace statefips = statefips[6] if mi(statefips)
assert ~mi(statefips)

merge m:1 statefips year using `mw', update
drop if _merge == 2 //Alaska, Hawaii, DC
assert ~mi(mw)
replace lmw = ln(mw) if mi(lmw)
assert ~mi(lmw)  
xtset pair_id_county year


frame copy mscz mscz_dl
frame change mscz_dl
forvalues i = 1985/1989 {
	insobs 1
	replace year = `i' if _n == _N
}
  
forvalues i = 2020/2022 {
	insobs 1
	replace year = `i' if _n == _N
}

fillin pair_id_czstate year
drop if mi(pair_id_czstate)
bys pair_id_czstate (year): replace statefips = statefips[6] if mi(statefips)
assert ~mi(statefips)

merge m:1 statefips year using `mw', update
drop if _merge == 2 //Alaska, Hawaii, DC
assert ~mi(mw)
replace lmw = ln(mw) if mi(lmw)
assert ~mi(lmw)  
xtset pair_id_czstate year


************************************************
* TWFE regressions
************************************************

foreach frame in twfe1 twfe1_00 twfe1_16 {
	frame `frame' {
		eststo `frame'_earn_uw: reghdfe learn75 lmw learn_other lpop, a(czstate year) cluster(statefips)
		eststo `frame'_emp_uw: reghdfe lemp75 lmw lemp_other lpop, a(czstate year) cluster(statefips)
		eststo `frame'_owe_uw: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop, a(czstate year) cluster(statefips)
		
		eststo `frame'_earn_w: reghdfe learn75 lmw learn_other lpop [aw=pop], a(czstate year) cluster(statefips)
		eststo `frame'_emp_w: reghdfe lemp75 lmw lemp_other lpop [aw=pop], a(czstate year) cluster(statefips)
		eststo `frame'_owe_w: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop [aw=pop], a(czstate year) cluster(statefips)

		estadd local fe " "
		estadd local czp_fe " "
		estadd local bcp_fe " "
		estadd local gap " "
		estadd local samples " "
		estadd local cz_state_sample "Y"
		estadd local mscz_sample " "
		estadd local bcp_sample " "	
	}
}



foreach frame in bcp bcp_00 bcp_16 {
	frame `frame' {
		eststo `frame'_earn_uw: reghdfe learn75 lmw learn_other lpop, a(pair_id_county pair_id_num##year) cluster(statefips)
		eststo `frame'_emp_uw: reghdfe lemp75 lmw lemp_other lpop, a(pair_id_county pair_id_num##year) cluster(statefips)
		eststo `frame'_owe_uw: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop, a(pair_id_county pair_id_num##year) cluster(statefips)
		
		eststo `frame'_earn_w: reghdfe learn75 lmw learn_other lpop [aw=pop], a(pair_id_county pair_id_num##year) cluster(statefips)
		eststo `frame'_emp_w: reghdfe lemp75 lmw lemp_other lpop [aw=pop], a(pair_id_county pair_id_num##year) cluster(statefips)
		eststo `frame'_owe_w: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop [aw=pop], a(pair_id_county pair_id_num##year) cluster(statefips)

		estadd local fe " "
		estadd local czp_fe " "
		estadd local bcp_fe "Y"
		estadd local gap " "
		estadd local samples " "
		estadd local cz_state_sample " "
		estadd local mscz_sample " "
		estadd local bcp_sample "Y"	
	}
}



foreach frame in mscz mscz_00 mscz_16 {
	frame `frame' {
		eststo `frame'_earn_uw: reghdfe learn75 lmw learn_other lpop, a(pair_id_czstate pair_id_num##year) cluster(statefips)
		eststo `frame'_emp_uw: reghdfe lemp75 lmw lemp_other lpop, a(pair_id_czstate pair_id_num##year) cluster(statefips)
		eststo `frame'_owe_uw: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop, a(pair_id_czstate pair_id_num##year) cluster(statefips)
		
		eststo `frame'_earn_w: reghdfe learn75 lmw learn_other lpop [aw=pop], a(pair_id_czstate pair_id_num##year) cluster(statefips)
		eststo `frame'_emp_w: reghdfe lemp75 lmw lemp_other lpop [aw=pop], a(pair_id_czstate pair_id_num##year) cluster(statefips)
		eststo `frame'_owe_w: ivreghdfe lemp75 (learn75 = lmw) learn_other lemp_other lpop [aw=pop], a(pair_id_czstate pair_id_num##year) cluster(statefips)
		
		estadd local fe " "
		estadd local czp_fe "Y"
		estadd local bcp_fe " "
		estadd local gap " "
		estadd local samples " "
		estadd local cz_state_sample " "
		estadd local mscz_sample "Y"
		estadd local bcp_sample " "	
	}
}

* Output twfe table
esttab twfe1_earn_uw mscz_earn_uw bcp_earn_uw twfe1_16_earn_uw mscz_16_earn_uw ///
bcp_16_earn_uw twfe1_00_earn_uw mscz_00_earn_uw bcp_00_earn_uw ///
using "$tables/Table 2.tex", replace booktabs b(3) se(3) nomtitles keep(lmw) ///
varlabels(lmw "Log wages") nonote noobs star(* 0.10 ** 0.05 *** 0.01) width(\hsize) ///
mgroups("1990-2019" "1990-2016" "2000-2016", pattern(1 0 0 1 0 0 1 0 0) span ///
prefix(\multicolumn{@span}{c}{) suffix(}) erepeat(\cmidrule(lr){@span})) ///
postfoot( ) ///
prehead(\begin{table}[h!]\centering ///
	\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}  ///
	\caption{MW effects on log average earnings and log employment in restaurants: TWFE estimates for different samples and specifications using QCEW data \label{tab:twfeqcew}}  ///
	\footnotesize ///
	\begin{tabular*}{\hsize}{@{\hskip\tabcolsep\extracolsep\fill}l*{9}{c}} \toprule) ///
posthead(\midrule \\ ///
	\multicolumn{10}{l}{\textbf{A. Unweighted}} \\ & & & \\)

esttab twfe1_emp_uw mscz_emp_uw bcp_emp_uw twfe1_16_emp_uw mscz_16_emp_uw ///
bcp_16_emp_uw twfe1_00_emp_uw mscz_00_emp_uw bcp_00_emp_uw /// 
using "$tables/Table 2.tex", append booktabs b(3) se(3) nomtitles keep(lmw) ///
varlabels(lmw "Log employment") nonote noobs nonum star(* 0.10 ** 0.05 *** 0.01) ///
width(\hsize) prehead( ) posthead( ) postfoot( )

esttab twfe1_owe_uw mscz_owe_uw bcp_owe_uw twfe1_16_owe_uw mscz_16_owe_uw ///
bcp_16_owe_uw twfe1_00_owe_uw mscz_00_owe_uw bcp_00_owe_uw /// 
using "$tables/table 2.tex", append booktabs b(3) se(3) nomtitles keep(learn75) ///
varlabels(learn75 "OWE") nonote noobs nonum star(* 0.10 ** 0.05 *** 0.01) width(\hsize) ///
prehead( ) posthead( ) postfoot( ) ///

esttab twfe1_earn_w mscz_earn_w bcp_earn_w twfe1_16_earn_w mscz_16_earn_w ///
bcp_16_earn_w twfe1_00_earn_w mscz_00_earn_w bcp_00_earn_w ///
using "$tables/Table 2.tex", append booktabs b(3) se(3) nomtitles keep(lmw) ///
varlabels(lmw "Log wages") nonote noobs nonum star(* 0.10 ** 0.05 *** 0.01) width(\hsize) ///
prehead(\multicolumn{10}{l}{\textbf{B. Weighted by population}} ///
\\ & & & \\) postfoot( ) posthead( )

esttab twfe1_emp_w mscz_emp_w bcp_emp_w twfe1_16_emp_w mscz_16_emp_w ///
bcp_16_emp_w twfe1_00_emp_w mscz_00_emp_w bcp_00_emp_w /// 
using "$tables/Table 2.tex", append booktabs b(3) se(3) nomtitles keep(lmw) ///
varlabels(lmw "Log employment") nonote noobs nonum star(* 0.10 ** 0.05 *** 0.01) ///
width(\hsize) prehead( ) posthead( ) postfoot( )

esttab twfe1_owe_w mscz_owe_w bcp_owe_w twfe1_16_owe_w mscz_16_owe_w ///
bcp_16_owe_w twfe1_00_owe_w mscz_00_owe_w bcp_00_owe_w ///
using "$tables/Table 2.tex", append booktabs b(3) se(3) nomtitles keep(learn75) ///
varlabels(learn75 "OWE") nonote noobs nonum star(* 0.10 ** 0.05 *** 0.01) width(\hsize) ///
stats(N gap fe czp_fe bcp_fe gap samples cz_state_sample mscz_sample bcp_sample, ///
labels("\textbf{N}" " " "\textbf{Fixed effects}" "MSCZ pair-year" "BC pair-year" ///
 " " "\textbf{Sample}" "All CZ-state" "MSCZ pairs" "Border county pairs") ///
fmt(%9.0fc)) /// 
prehead( ) posthead( ) ///
postfoot( ///
\hline\hline ///
\end{tabular*} ///
{\centering ///
\caption*{\begin{footnotesize} ///
\textit{Notes:} Estimated using specification described in equation \ref{eq:TWFE_logMW} and Quarterly Census of Employment and Wages data. Columns (1)-(3) use the years 1990-2019, (4)-(6) use 1990-2016, (7)-(9) use 2000-2016. Columns (1), (4), and (7) use the sample of all commuting-zones-by-state, and include CZ-state and year fixed effects. Columns (2), (5), and (8) use multi-state commuting-zone-by-state pairs sample and include fixed effects for CZ-state (repeated if a CZ-state is part of multiple pairs) and pair-year fixed effects. Columns (3), (6), and (9) similarly use the sample of border county pairs and include fixed effects for border counties (each time a county is in a pair) and pair-year fixed effects. Outcomes include log weekly earnings, log employment, and own wage elasticity (OWE). OWE is the IV estimate from regressing log employment on log earnings instrumented by log minimum wage. OWE estimates are only reported when there is a strong first-stage effect. Estimates in Panel A are unweighted, while estimates in Panel B are weighted using the working age population. All specifications include log earnings/log employment outside of the restaurant sector, and working-age population as control variables. Standard errors, reported in parentheses, are clustered at the state level. Stars indicate statistical significance as follows: sym{*} \(p<0.10\), \sym{**} \(p<0.05\), \sym{***} \(p<0.01\) ///
\end{footnotesize}}} ///
\end{table} ///
) ///
substitute("\_" "_")



************************************************
* Rolling regressions
************************************************

* max samples (till 2016 and 2019) with different start years (1990-2005)
frame change mscz

foreach j in 2016 2019 {
	forvalues i = 1990/2005 {
		
		reghdfe lemp75 lmw lemp_other lpop if inrange(year,`i',`j') [aw=pop], ///
		a(pair_id_czstate pair_id_num##year) cluster(statefips)
		
		if `i' == 1990 {
			matrix w_estimates`j' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
		}
		
		else {
			matrix w_estimates`j'_`i' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
			matrix w_estimates`j' = (w_estimates`j' \ w_estimates`j'_`i')
		}
		
		reghdfe lemp75 lmw lemp_other lpop if inrange(year,`i',`j') [aw=lpop], ///
		a(pair_id_czstate pair_id_num##year) cluster(statefips)
		
		if `i' == 1990 {
			matrix w_lpop_estimates`j' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
		}
		
		else {
			matrix w_lpop_estimates`j'_`i' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
			matrix w_lpop_estimates`j' = (w_lpop_estimates`j' \ w_lpop_estimates`j'_`i')
		}
		
		
		reghdfe lemp75 lmw lemp_other lpop if inrange(year,`i',`j'), ///
		a(pair_id_czstate pair_id_num##year) cluster(statefips)
		
		if `i' == 1990 {
			matrix uw_estimates`j' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
		}
		
		else {
			matrix uw_estimates`j'_`i' = (r(table)[1,1], r(table)[5,1],r(table)[6,1])
			matrix uw_estimates`j' = (uw_estimates`j' \ uw_estimates`j'_`i')
		}
		
	}

foreach mat in w_estimates`j' w_lpop_estimates`j' uw_estimates`j' {
	matrix colnames `mat' = b_qcew ll_qcew ul_qcew
}

	
}

foreach i in w uw w_lpop {
	matrix estimates_`i' = (`i'_estimates2016 \ `i'_estimates2019)
	frame create estimates_`i'
	frame change estimates_`i'
	svmat estimates_`i', names(col)
	gen year = _n + 1989 if _n <= 16
	replace year = _n + 1973 if _n > 16
	foreach var in b_qcew ll_qcew ul_qcew {
		gen `var'_2019 = `var' if _n > 16
		replace `var' = . if _n > 16
		rename `var' `var'_2016
	}
	append using "$out/intermediate_results/rreg_cbp_estimates_`i'.dta"
	replace year = year + 0.2 if ~mi(b_qcew_2019)
	replace year = year - 0.2 if ~mi(b_cbp)
}


* Outputting results

frame change estimates_w

#delimit ;
twoway 
	(scatter b_cbp year, color(red*0.5) msymbol(s)) 
	(rcap ll_cbp ul_cbp year, lcolor(red*0.25) lwidth(medthin)) 
	(scatter b_qcew_2016 year, color(blue*0.5) msymbol(o)) 
	(rcap ll_qcew_2016 ul_qcew_2016 year, lcolor(blue*0.25) lwidth(medthin)) 
	(scatter b_qcew_2019 year, color(green*0.5) msymbol(oh)) 
	(rcap ll_qcew_2019 ul_qcew_2019 year, lcolor(green*0.25) lwidth(medthin)),
	xscale(range(1989.6 2005.4))
	xlab(1990(5)2005.4, labsize(small)) 
	yscale(range(-0.52 0.2))
	ylab(-0.5(0.1)0.2, labsize(small))
	yline(0, lcolor(gs12) lpat(solid)) 
	xtitle(Sample start year, size(small)) 
	ytitle(MW elasticity, size(small)) 
	legend(order(1 3 5) label(1 "CBP up to 2016") label(3 "QCEW up to 2016") label(5 "QCEW up to 2019")  pos(6) rows(1)) 
	name(weighted, replace) 
	subtitle("B. Weighted by population", size(small))
	;
#delimit cr

frame change estimates_uw

#delimit ;
twoway 
	(scatter b_cbp year, color(red*0.5) msymbol(s)) 
	(rcap ll_cbp ul_cbp year, lcolor(red*0.25) lwidth(medthin)) 
	(scatter b_qcew_2016 year, color(blue*0.5) msymbol(o)) 
	(rcap ll_qcew_2016 ul_qcew_2016 year, lcolor(blue*0.25) lwidth(medthin)) 
	(scatter b_qcew_2019 year, color(green*0.5) msymbol(oh)) 
	(rcap ll_qcew_2019 ul_qcew_2019 year, lcolor(green*0.25) lwidth(medthin)),
	xscale(range(1989.6 2005.4))
	xlab(1990(5)2005.4, labsize(small)) 
	yscale(range(-0.52 0.2))
	ylab(-0.5(0.1)0.2, labsize(small))
	yline(0, lcolor(gs12) lpat(solid)) 
	xtitle(Sample start year, size(small)) 
	ytitle(MW elasticity, size(small)) 
	legend(off) 
	name(unweighted, replace) 
	title(" ")
	subtitle("A. Unweighted", size(small))
	;
#delimit cr


grc1leg2 unweighted weighted, legendfrom(weighted) rows(2) ysize(7) xsize(7) ///
	imargin(10 10 0 0)

graph export "$figures/Figure 1.pdf", replace


frame change estimates_w_lpop

#delimit ;
twoway 
	(scatter b_cbp year, color(red*0.5) msymbol(s)) 
	(rcap ll_cbp ul_cbp year, lcolor(red*0.25) lwidth(medthin)) 
	(scatter b_qcew_2016 year, color(blue*0.5) msymbol(o)) 
	(rcap ll_qcew_2016 ul_qcew_2016 year, lcolor(blue*0.25) lwidth(medthin)) 
	(scatter b_qcew_2019 year, color(green*0.5) msymbol(oh)) 
	(rcap ll_qcew_2019 ul_qcew_2019 year, lcolor(green*0.25) lwidth(medthin)),
	xscale(range(1989.6 2005.4))
	xlab(1990(5)2005.4, labsize(small)) 
	yscale(range(-0.52 0.2))
	ylab(-0.5(0.1)0.2, labsize(small))
	yline(0, lcolor(gs12) lpat(solid)) 
	xtitle(Sample start year, size(vsmall)) 
	ytitle(MW elasticity, size(vsmall))
	name(unweighted, replace) 
	title(" ")
	subtitle("Weighted by log population", size(small))
	legend(order(1 3 5) label(1 "CBP up to 2016") label(3 "QCEW up to 2016") label(5 "QCEW up to 2019")  pos(6) rows(1)) 
	;
#delimit cr


graph export "$out/figures/Figure A3.pdf", replace



************************************************
* TWFE DL figures
************************************************


frame change mscz_dl

reghdfe lemp75 L(1/5).lmw lmw F(1/3).lmw lemp_other lpop if inrange(year,1990,2016), a(pair_id_czstate pair_id_num##year) cluster(statefips)

nlcom (t_4: -(_b[F1.lmw] + _b[F2.lmw] + _b[F3.lmw])) ///
(t_3: -(_b[F1.lmw] + _b[F2.lmw])) (t_2: -_b[F1.lmw]) (t_1: 0) ///
(t0: _b[lmw]) (t1: _b[lmw] + _b[L1.lmw]) ///
(t2: _b[lmw] + _b[L1.lmw] + _b[L2.lmw]) ///
(t3: _b[lmw] + _b[L1.lmw] + _b[L2.lmw] + _b[L3.lmw]) ///
(t4: _b[lmw] + _b[L1.lmw] + _b[L2.lmw] + _b[L3.lmw] + _b[L4.lmw]) ///
(t5: _b[lmw] + _b[L1.lmw] + _b[L2.lmw] + _b[L3.lmw] + _b[L4.lmw] + _b[L5.lmw])

matrix step1 = r(table)[1,1..10]
matrix step2 = r(table)[5..6,1..10]
matrix lemp = (step1 \ step2)
matrix lemp = lemp'

frame create mscz_dl_lemp
frame change mscz_dl_lemp
svmat lemp, names(col)
foreach v in ll ul {
	replace `v'=0 if mi(`v')
}
gen eventtime = _n - 5


frame mscz_dl_lemp: twoway (line b eventtime, lcolor(blue*0.5) lpattern(solid)) ///
(rcap ll ul eventtime, lcolor(blue*0.5) lpattern(dash)), ///
xtitle(Years around MW increase event) xlabel(-4 "{&le} -4" -3(1)4 5 "{&ge} 5") ///
ytitle("Cumulative response of log employment") legend(off) ///
yline(0) ylabel(-0.8(0.2)0.4) title(A. Sample - QCEW) ///
saving("$figures/mscz_qcew_dl", replace)

gr combine "$figures/mscz_qcew_dl.gph" "$figures/mscz_cbp_dl.gph", rows(1) xsize(7) ysize(3)
graph export "$figures/Figure A4.pdf", replace

